# load iris example dataset -------------------------------------------------------
data(iris)


# inspect dataset ---------------------------------------------------------
iris
iris$Sepal.Length

iris[1:3]    #returns first three columns
iris[1:3,]   #returns first three rows
iris[25:27,2:4]  #returns first three columns for 5th row

head(iris)
head(iris[iris$Species=="virginica",])
str(iris)



# descriptive statistics --------------------------------------------------
summary(iris)
summary(iris[iris$Species=="versicolor",])
mean(iris$Sepal.Length)
sd(iris$Sepal.Length)
median(iris$Petal.Width)
var(iris$Sepal.Length[iris$Species=="virginica"])


attach(iris)
mean(Sepal.Length)

hist(Sepal.Length)

library(ggplot2)
ggplot(data=iris, aes(x=Sepal.Length, fill=Species)) +
  geom_histogram(binwidth=0.2, alpha=0.5, position="identity")

detach(iris)
detach()



#install.packages("pastecs")
library(pastecs)
stat.desc(iris)
stat.desc(iris)[11,2]


#install.packages("psych")
library(psych)

attach(iris)
describeBy(iris, group=Species)
describeBy(iris, group=Species)[3]
?describeBy


#effect size 
?cohen.d

cohen.d(iris$Sepal.Length, iris$Species==c("virginica","versicolor")) #finds the standardized mean difference (Cohen's d)
res_d<-cohen.d(iris$Sepal.Length, iris$Species==c("virginica","versicolor")) #finds the standardized mean difference (Cohen's d)
res_d
error.dots(res_d)




# correlation -------------------------------------------------------------
Sepal.length.setosa<-iris$Sepal.Length[iris$Species=="setosa"]
Sepal.width.setosa<-iris$Sepal.Width[iris$Species=="setosa"]

plot(Sepal.length.setosa,Sepal.width.setosa , xlab="Sepal Length",ylab="Sepal width", main="Iris setosa")

#assumptions - normality
hist(Sepal.length.setosa)
hist(Sepal.width.setosa)

shapiro.test(Sepal.length.setosa)
shapiro.test(Sepal.width.setosa)

#correlation tests
cor(Sepal.length.setosa,Sepal.width.setosa)

cor.test(Sepal.length.setosa,Sepal.width.setosa)
?cor.test

library(psych)
pairs.panels(iris, method="pearson")




# tests for differences in means ------------------------------------------

# t-tests

Sepal.length.setosa<-iris$Sepal.Length[iris$Species=="setosa"]
Sepal.length.versicolor<-iris$Sepal.Length[iris$Species=="versicolor"]

boxplot(Sepal.length.setosa, Sepal.length.versicolor, names=c("setosa","versicolor"), ylab="Sepal Length")


#assumptions normal distribution and variance of the two groups is equal (homoscedasticity)
shapiro.test(Sepal.length.setosa)
shapiro.test(Sepal.length.versicolor)
var.test(Sepal.Length~Species, data=iris[iris$Species==c("virginica","versicolor"),])


#one sample t-test
t.test(Sepal.length.setosa, mu=5)
t.test(Sepal.length.setosa, mu=3)


#two-sample t-test (assumes unequal variance)
t.test(Sepal.length.setosa, Sepal.length.versicolor)
?t.test

# Analysis of Variance - ANOVA
boxplot(Sepal.Length~Species, data=iris)
fit<-aov(Sepal.Length~Species, data=iris)
summary(fit)

#Tukey Honestly Significant Differences
TukeyHSD(fit)



# linear models -----------------------------------------------------------
data(mtcars)

head(mtcars)
str(mtcars)

#install.packages("car")
library(car)

#simple linear model
mod1<-lm(mpg~hp, data=mtcars)
mod1
plot(mtcars$mpg~mtcars$hp)
abline(mod1)
summary(mod1)
coefficients(mod1)




# bonus material ----------------------------------------------------------


#multiple regression
mod2<-lm(mpg~disp+hp+wt+drat, data=mtcars)
summary(mod2)

coefficients(mod2)

#diagnostic plots and Outlier test
par(mfrow=c(2,2))
plot(mod1)
outlierTest(mod1)

#comparing model fits
mod3<-lm(mpg~hp+wt, data=mtcars)
summary(mod3)
AIC(mod2, mod3)



install.packages("MASS")
library(MASS)
step <- stepAIC(mod2, direction="both")
step$anova # display results




# estimate power of a test - important for experimental design ------------


library(psych)
#cohen's d = standardized mean difference
Sepal.cohen<-cohen.d(iris$Sepal.Length, iris$Species==c("virginica","versicolor"))$cohen.d[2]
Sepal.sd<-sd(c(Sepal.length.setosa, Sepal.length.versicolor))

power.t.test(n = 50, d = NULL, sd=Sepal.sd, sig.level = 0.05, power = 0.9, type = "two.sample")
?power.t.test


# checking the power test
mu1<-5
mu2<-6
s<-1.2
n<-12
X1<-rnorm(n, mu1, s)
X2<-rnorm(n, mu2, s)

power.t.test(n,delta=mu2-mu1, sd=s)

Pv<-vector()
for(i in 1:10000){
  X1<-rnorm(n, mu1, s)
  X2<-rnorm(n, mu2, s)
  my.test<-t.test(X1,X2, var.equal = T)
  Pv[i]<-my.test$p.value}
par(mfrow=c(1,1))
hist(Pv)
sum(Pv<0.05)